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ABSTRACT 


A  fifteen-month  multi-tasked  research  project  was  pursued  by  the 


present  investigators  to  study  complex  viscous  flows  under  AFOSR  sponsorship 
between  July  1985  and  September  1986.  The  major  objective  of  this  study  was 


to  acquire  improved  understanding  of  viscous  flows  and  to  develop  basic 


computational  methods  for  efficient  determination  of  2-D/3-D  subsonic  and 


incompressible  flows.  Two  major  analyses  were  pursued.  These  include  the 
Interacting  Parabolised  Navier-Stokes  (IPNS)  analysis  for  steady  flows  and 


the  full  Navier-Stokes  (NS)  analysis  for  direct  simulation  of  unsteady 


flows.  The  IPNS  analysis  developed  employs  no  ad  hoc  artificial  dissipation 


and,  in  spite  of  being  a  density-based  formulation,  performs  well  even  for 


very  low  Mach  numbers.  The  applications  considered  include  2-D  cascades  and 


channels  of  simple  geometry.  The  flow  solutions  are  well  behaved  in  the 


presence  of  sharp  leading  edges  and  trailing  edges,  as  well  as  in  the 

i 

presence  of  reversed  flow.  The  complete  unsteady  NS  analysis  is  based  on 


direct- solution  methodology  and  has  proved  to  be  very  robust  and  efficient. 


It  has  been  applied  to  analyze  high- incidence  aerodynamics  of  symmetric  as 


well  as  cambered  Joukowski  airfoils,  and  has  yielded  very  interesting 


results  regarding  multiple  incommensurate  frequencies  in  the  large-time 


behavior  of  the  flow  and,  hence,  regarding  the  theory  of  strange  attractors 

' 

i 

I  and  chaos.  The  corresponding  3-D  analysis  appears  wall  on  its  way  and  looks 


very  promising.  All  of  this  research  relies  heavily  on  usage  of  the 


supercomputers  at  the  NASA  Research  Centers,  in  addition  to  all  the  various 


computing  and  graphics  facilities  at  the  University  of  Cincinnati.  Both  the 


2-D  analyses  developed  are  reasonably  mature  now  and  should  be  useful  for 


examining  realistic  flow  phenomena. 
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SECTION  1 


OBJECTIVES 

The  development  of  computational  fluid  dynamics  (CFD)  analyses  for 
complex  viscous  interacting  flows  was  pursued  by  the  present  investigators 
during  July  1985  -  September  1986.  The  major  thrust  of  the  work  was  in  the 
development  of  two  analyses,  which  are  briefly  outlined  as  under: 

Viscous  Interacting  Analysis:  Internal  Flows 

Subsonic  flows  are  essentially  elliptic  in  character.  However,  in  a 
large  class  of  internal  flow  applications,  streamwise  diffusion  is 
negligible,  although  upstream  Influences  through  pressire  interactions  are 
significant.  These  situations  can  be  analyzed  and  computed  more  efficiently 
using  the  Interacting  Parabolized  Navier-Stokes  (IPNS)  model  rather  than  the 
full  Navier-Stokes  equations,  particularly  for  steady  flow.  The  IPNS  model 
and  the  corresponding  numerical  method  were  studied,  with  the  objective  of 
developing  an  efficient  procedure  of  desirable  versatility  for  use  in 
studying  flows  in  diffusers,  cascades,  etc. 

Unsteady  Navier-Stokes  Analysis:  Internal  and  External  Flows 

Unsteady  separation  and  vorticity  interactions  in  self-excited  and 
forced  unsteady  flows  are  important  fundamental  phenomena  which  occur  in 
internal  as  well  as  external  flows.  Unsteady  separation  in  large-vortex 
dominated  flows  is  particularly  important  in  external  flows.  The  underlying 
physical  mechanisms  were  examined  by  analyzing  flows  exhibiting  these 
phenomena,  with  the  help  of  the  complete  unsteady  Navier-Stokes  equations 
solved  using  a  direct  (fully  implicit)  numerical  technique. 


Additional  items  such  as  grid  generation  and  development  of  solution 
algorithms  were  also  studied  to  further  aid  the  above  two  analyses.  These 
items  contribute  significantly  to  basic  research  even  by  themselves. 
Moreover,  these  essential  elements  were  needed  for  successfully  developing 
the  two  main  analyses. 


SECTION  2 


DESCRIPTION  OF  SIGNIFICANT  ACCOMPLISHMENTS 
All  of  the  areas  of  research  proposed  for  study  were  pursued  during 
this  15-month  grant  period.  The  progress  achieved  in  each  of  these  areas  is 
briefly  described  in  this  section. 

2. 1  Composite  Generation  of  Multi-Block  Grids  for  Subsonic  Cascades 

Examination  of  the  basic  grids  for  cascades  shows  that  a  C-grid  can 
adequately  resolve  the  boundary  layers  on  the  blades  as  well  as  the  wakes 
downstream  and  is  satisfactory  everywhere  except  in  the  region  upstream  of 
the  blades.  In  this  region,  the  grid  density  decreases  rapidly  with 
distance  upstream  of  the  front  stagnation  point.  On  the  other  hand,  an 
H-grid  is  generally  satisfactory  in  this  region  and  also  facilitates 
implementation  of  the  repeating  boundary  condition  for  cascade  flows.  For 
cascades  with  minimal  stagger  and  blades  with  thin  leading  edges,  the  H-grid 
is  not  highly  skewed.  But  when  the  blades  possess  thick  rounded  leading 
edges  typical  of  turbine  cascades,  the  H-grid  becomes  highly  skewed  in  the 
leading-edge  regions  (Fig.  1).  Earlier,  K.  Ghia  and  U.  Ghia  [1982]  had 
suggested  an  approach  for  improving  a  basic  C-grid  for  use  with  cascades. 
This  consisted  of  deleting  a  segment  of  the  C-grid  in  the  upstream  region 
and  patching  a  segment  of  an  H-grid  onto  the  remainder  of  the  C-grid,  as 
shown  in  Fig.  2.  The  resulting  coordinates  appear  satisfactory  everywhere. 
It  should  be  noted  that  this  hybrid  coordinate  system  contains  two  five¬ 
sided  cells  which  require  special  consideration. 

The  concept  of  C-H  hybrid  grids  for  cascades  was  pursued  further  by 
U.  Ghia  et  al .  [1983].  The  main  result  of  that  effort  is  shown  in  Fig.  3. 
Here,  a  C-grid  encompasses  a  narrow  region  in  the  vicinity  of  the  lower 


blade  and  approximately  includes  the  wake  region.  Similarly,  a  C-grid  also 


covers  the  corresponding  repeating  region  around  the  upper  blade.  The 
remainder  of  the  flow  domain  between  the  blades  is  then  occupied  by  a 
sequence  of  three  H-grid  blocks.  This  led  to  a  five-block  grid  system  for 
the  cascade.  The  grid  in  each  of  the  five  blocks  was  generated  separately 
using  numerical  solution  of  elliptic  partial  differential  equations 
governing  the  coordinate  transformation.  Continuity  across  the  block 
interfaces  was  maintained  by  a  single  visit  to  each  block,  performed  in  a 
specific  sequence.  First,  the  C-grid  portions  of  the  cascade  coordinates 
were  generated,  using  appropriate  forcing  functions  to  provide  resolution  of 
the  blade  boundary  layers.  The  resultant  clustering  near  the  outer  boundary 
of  the  C-grid  blocks  was  employed  to  determine  the  forcing  functions  for  the 
grids  in  the  adjoining  H-grid  blocks,  leading  to  an  overall  grid  system  that 
is  continuous  everywhere.  This  procedure  also  circumvented  the  difficulties 
to  be  expected  with  the  two  five-sided  cell3  appearing  in  this  hybrid  grid. 
However,  the  procedure  required  knowledge  of  the  grid-point  distribution  at 
all  block  interfaces;  this  constituted,  a  major  difficulty,  especially  so  far 
as  the  solution  for  the  flow  itself  was  concerned. 

Work  performed  under  the  present  grant  ha3  resulted  in  the  development 
of  a  composite  solution  procedure  which  does  not  employ,  or  require,  any 
information  along  block  interfaces.  The  boundaries  of  the  computational 
region  correspond  only  to  actual  boundaries  in  the  physical  plane.  All 
information  along  physical-block  interfaces  is  obtained  as  part  of  the 
evolving  overall  solution.  Therefore,  no  special  iterative  measures  are 
required  at  these  interfaces.  The  first  step  towards  developing  this 
procedure  was  the  appropriate  representation  of  the  physical  region  in  the 
computational  domain.  The  earlier  3"D  computational- domain  representation 


Uk 


(Fig.  3)  by  U.  Ghia  et  al.  [1983]  of  the  2-D  hybrid  grid  was  very  useful  but 
was  computationally  inefficient.  The  present  representation  of  the  2-D 
physical  problem  retains  a  2-D  computational  region.  Only  the  true 
boundaries  of  the  physical  region  appear  as  boundaries  in  the  computational 
plane;  no  block  interfaces  are  exposed  as  computational  boundaries.  The 
resultant  multi-rectangular  computational-region  representation  also  paves 
the  way  for  appropriate  treatment  of  the  five-sided  cells  appearing  in  the 
hybrid  grid.  Finally,  this  representation  is  very  well  suited  for  the 
numerical  solution  of  the  coordinate-generation  equations,  as  well  as  the 
flow  equations,  even  by  implicit  methods  such  as  ADI.  The  various  blocks 
used  in  the  hybrid  grid  for  a  staggered  cascade  are  shown  in  Fig.  4  which 
shows  mainly  the  boundaries  of  these  blocks.  The  corresponding 
computational  region  is  shown  in  Fig.  5.  With  this  arrangement,  the  grid  is 
developed  simultaneously  in  all  the  blocks  without  any  special  consideration 
of  block  interfaces.  This  is  a  unique  feature  of  the  present  multi-block 
domain-decomposition  procedure.  An  invited  paper  based  on  some  preliminary 
results  from  this  work  was  presented  by  U.  Ghia  et  al.  [1985];  a  written 
version  of  the  paper  was  also  prepared  later  by  U.  Ghia  et  al.  [1986]. 

2.2  Steady  Flows  with  Strong  Upstream  Interactions 

In  this  effort,  a  reduced  form  of  the  governing  equations  has  been 
developed  which  can  capture  much  of  the  physics,  while  requiring  less 
computer  resources  than  the  full  Navier-Stokes  equations.  It  belongs  to  the 
category  of  semi-elliptic  analyses,  one  form  of  which  was  proposed  and 
employed  earlier  by  U.  Ghia  et  al.  [1981].  The  formulation  holds  greatest 
promise  for  high-Re  steady  flows  with  a  predominant  flow  direction.  It 
should  be  mentioned,  however,  that  the  procedure  can  be  extended  easily  to 


include  consideration  of  unsteady  flows.  But  the  most  unique  feature  of  the 
nethod  is  its  ability  to  conpute  low-Mach  nunber  flows.  Although  it  is  a 
density-based  method,  it  is  not  plagued  by  computational  difficulties  for 
M  -*  0,  as  are  the  other  density-based  methods  available.  Also,  no  use  is 

made  of  any  externally  added  artificial  viscosity.  The  method  is  applicable 
to  3-D  flows  as  well. 

The  governing  differential  equations  used  are  derived  from  the  Navi er- 
Stokes  equations,  together  with  the  continuity  equation  and  the  energy 
equation,  for  steady  flow  of  a  compressible  fluid.  With  Cartesian 
decomposition  of  the  vector  quantities,  these  equations  are  expressed,  in 
terms  of  a  general  coordinate  system  (£,n) ,  in  the  strong-conservation  law 
(SCL)  form.  Expressing  the  transformed  governing  differential  equations  in 
the  SCL  form  requires  that  the  coordinate  transformation  metrics  satisfy  a 
set  of  geometric  conservation  relations.  In  an  analytical  formulation, 
these  are  identically  satisfied.  In  a  discretized  formulation,  satisfaction 
of  these  relations  is  ensured  by  the  use  of  appropriate  differencing  for  the 
metrics. 

The  semi-elliptic  formulation  developed  in  the  present  work  is  obtained 
by  invoking  the  approximation  that,  for  flows  with  a  predominant  flow 
direction,  streamwlse  diffusion  is  negligible  relative  to  normal  diffusion. 
Nevertheless,  a  large  class  of  flows,  for  which  streamwlse  diffusion  may 
well  be  negligible,  are  significantly  Influenced  by  upstream  Interactions, 
and  may  not  be  adequately  represented  by  mathematically  parabolic  equations. 
In  the  present  formulation,  upstream  interaction  is  provided  for  via 
appropriate  treatment  of  the  streamwlse  pressure  gradient  term.  Hence,  the 
formulation  1s  also  termed  an  interacting  parabolized  Navler-Stokes  (IPNS) 


model.  The  viscous- lnviscid  interaction  is  Included  by  composite 


consideration  of  the  viscous  and  inviscid  flow  regions  through  the  use  of 
the  PNS  equations;  upstream  interactions  are  propagated  and  included  through 
the  pressure  field. 

For  the  subsonic  flows  considered,  the  total  pressure,  total 

temperature  and  streamline  slope  are  prescribed  at  Inflow  and  the  static 

pressire  is  prescribed  at  outflow.  At  the  wall-wall  boundaries,  the 

conditions  of  zero  slip  and  zero  suction/ injection,  together  with  a 

specified  temperature  condition,  provide  a  total  of  six  boundary  conditions. 

The  one  additional  boundary  condition  needed  is  provided  by  an  approximate 

form  of  the  normal  momentum  equation  obtained  by  neglecting  the  viscous 

terms  in  that  equation.  This  last  condition  is  applied  in  the  cells 

adjacent  to  the  wall,  rather  than  at  the  wall  itself,  thus  avoiding  the  need 

for  any  one-sided  differences  for  the  normal  derivatives.  The  wake-wake 

boundaries  are  the  periodic  boundaries  occurring  in  cascade  flows.  The 

periodicity  condition  requires  that  the  corresponding  values  of  all  four 

flow  variables,  and  the  rr derivatives,  u  ,  v  and  T  ,  of  the  velocities  and 

n  n  n 

temperature  which  are  governed  by  second-order  differential  equations,  be 
the  same  at  corresponding  periodic  points  along  the  wake  boundaries.  It  is 
Important  to  mention  that,  in  terms  of  the  conserved  variables  (p,  pu,  pv, 

pet)  comprising  the  solution  vector  Q,  the  repeating  condition  on  the 

n-derivatlves  must  be  satisfied  for  all  four  elements  of  Q  . 

n 

The  analysis  was  employed  to  determine  the  flow  in  several  channel  and 
cascade  configurations.  Some  of  the  results  obtained  are  presented  here  in 
terms  of  the  two  most  sensitive  quantities,  namely,  the  static  pressure  and 
wall  shear.  Unless  otherwise  stated,  the  Mach  number  M  is  approximately 


N 


Figure  6  shows  the  distribution  of  the  surface  pressure 

p  (-  p  -p  )  and  wall-shear  parameter,  t  (-  3u  /3n),  obtained  for  a 
D  winlet  w  w  s 

constricted  channel  with  t _ /h  -  0.2  and  Re  «  1500.  Of  particular 

max 

significance  is  the  presence  of  a  finite  region  of  small  reversed  flow 
indicated  by  negative  tw  downstream  of  the  constriction  and  the  completely 

non-singular  behavior  of  the  present  IPNS  solution  for  this  configuration 
with  upstream  interaction. 

The  performance  of  the  IPNS  model  in  the  presence  of  sharp  leading  and 
trailing  edges  was  evaluated  by  application  of  the  model  to  a  cascade  of 
finite  flat  plates.  Figure  7  shows  the  distribution  of  the  pressure  pfe  and 

shear  parameter  xw  for  Re  ranging  from  1500  to  15,000.  Through  all  the 

highly  nonlinear  behavior  of  the  flow  variables,  including  that  due  to 
abrupt  change  in  the  boundary  conditions  across  the  leading  edge  (LE)  and 
the  trailing  edge  (TE),  the  solution  of  the  IPNS  model  is  quite  regular  as 
upstream  interactions  are  appropriately  included  in  it. 

The  distributions  of  p^  and  xy  are  shown  in  Fig.  8  for  parabolic-arc 

airfoils  for  various  thickness  ratios  (Fig.  8a)  and  for  various  Reynolds 
numbers  (Fig.  8b).  The  effect  of  varying  M  up  to  nearly  0.5  has  also  been 

examined.  It  is  found  to  be  minimal  on  the  shear  parameter  tw,  but  becomes 

evident  in  the  pressure  distribution  at  the  highest  value  of  M-  considered. 

Figure  8c  shows  the  static  pressure  contours  for  the  case  with  M-  -  0.49,  Re 


-  15,000.  The  high-pressure  region  localized  near  the  LE  is  evident  from 
the  concentration  of  the  contours  in  this  region,  as  is  the  more  widespread 
low-pressure  regions  downstream  of  the  position  of  maximun  thickness.  From 
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the  values  of  the  pressure  along  these  contours,  it  is  clear  that  the 
pressure  varies  rather  minimally  over  the  airfoil  surface.  Nevertheless, 
these  variations  have  been  very  accurately  computed  in  order  to  produce  the 
contours  which  are  well  behaved  and  conform  to  the  physics  of  this  flow, 
especially  in  the  LE  and  TE  regions. 

Application  to  these  various  flow  configurations  served  to  demonstrate 
that  the  technique  is  viable  for  flows  with  strong  interactions  occurring 
due  to  boundary^layer  separation  or  the  presence  of  sharp)  leading/trailing 
edges.  It  is  important  to  recall  that  the  procedure  developed  uses  no 
externally  added  artificial  viscosity  and  is  capable  of  producing 
satisfactory  solutions  for  compressible  (subsonic)  viscous  flows,  with  no 
modification  needed  for  analyzing  nearly  incompressible  flow  as  well. 

A  paper  based  on  these  results  was  presented  by  U.  Ghia  et  al.  [1985]; 
a  written  version  of  the  paper  was  also  prepared  later  by  U.  Ghia  et  al. 
[1986]. 

2.3  Unsteady  2-D  Navier-Stokes  Analysis 

The  simulation  of  2-D  flows  with  self- sustained  unsteadiness  has  been 
continued  using  the  direct  solution  of  the  unsteady  2-D  Navier-Stokes 
equations.  In  addition  to  determining  the  basic  flow  solution,  the  time- 
dependent  aerodynamic  lift,  drag  and  moment  coefficients  were  also  obtained 
for  flow  past  airfoils  at  high  Incidence.  Effort  was  also  made  to 
understand  the  observed  quasi  periodicity  and  bring  forth  any  possible 
similarity  with  strange  attractors. 


Symmetric  Joukowskl  Airfoil 

A  12  percent  thick  symmetric  Joukowskl  airfoil  is  used  in  this  study  as 
it  has  two  especially  attractive  features,  (i)  The  Joukowskl  airfoil  can  be 
accurately  represented  using  conformed  transformations;  the  details  of  these 
and  the  clustering  transformations  used,  as  well  as  the  total  analysis,  were 
given  by  K.  Ghia  et  al.  [1985a].  (ii)  The  presence  of  a  sharp  TE  leads  to  a 
much  stronger  interacting  region  and,  hence,  truly  tests  the  analysis 
developed.  This  unsteady  Navi er-S tokes  analysis  and  the  corresponding 
numerical  method  are  used  to  study  three  flow  configurations  in  detail.  All 
of  these  configurations  have  the  same  Re  -  1 ,000  but  the  value  of  the  free- 
stream  incidence  angle  varies  such  that  of  •  15°.  30°  and  53°. 

respectively.  For  the  high  angle-of-attack  case  with  af«  53#.  the  Joukowskl 

airfoil  appears,  to  the  oncoming  stream,  as  an  apparent  bluff  body.  The 
massively  separated  vortex-dominated  flow  in  the  post-stall  regime  for  this 
configuration  of  the  Joukowskl  airfoil  is  exceedingly  complex  and,  from  the 
results  available  so  far,  it  is  feasible  to  conjecture  two  hypotheses;  see 
Fig.  11.  One  possibility  is  that  the  solution  has  still  not  asymptoted  to 
an  exact  limit  cycle  but  may  do  so  subsequently.  The  second  possibility, 
based  on  the  results  for  af  -  30°,  is  that  the  solution  may  asymptote  to  a 

quasiperiodic  state,  with  anywhere  from  3  to  8  incommensurate  frequencies. 
The  second  hypothesis  is  more  likely  to  prevail. 

Limit  Cycle  Analysis 

Originally,  the  aerodynamic  coefficients  were  computed  only  at 
intervals  of  0.1  characteristic  time.  This  was  quite  satisfactory  for 
qualitative  assessment  of  the  flow  evolution  but  too  coarse  for  its  detailed 


analysis.  Subsequently,  the  computer  program  has  been  modified  to  provide 


the  coefficients  of  lift  C,  ,  drag  C_  and  moment  Cu  (nose-up  positive  about 

L  D  M 

the  quarter-chord  point)  at  every  At  increment.  The  configuration  with 
af  -  15°,  which  requires  minimum  time  to  reach  the  time-asymptotic  limit 

cycle  solution,  has  been  completely  recomputed,  using  a  slightly  improved 
grid  near  the  TE.  The  flow  configuration  with  -  30°  has  been  recomputed 

between  t  -  M5  and  t  -  58,  whereas,  due  to  limitation  of  availability  of  CPU 
time  on  the  host  computer  system,  the  configuration  with  af  -  53°  is 

currently  available  with  C,  ,  C_  and  Cu  computations  at  time  intervals  of  0.1 

L  D  M 

only.  Figures  9a,b,c  show  C,  ,  Cn  and  Cu  corresponding  to  cu«15°.  As  seen 

L  U  M  r 

in  this  figure,  C.  rises  initially  but  drops  very  sharply  during  the 
L 

transient  phase  and  asymptotes  to  a  near-limit-cycle  solution  corresponding 
to  the  dominant  frequency  for  the  shedding  of  vortices  from  the  TE.  The 
L2-norms  of  the  entire  vorticity  and  stream-function  fields  were  carefully 

examined  to  ensure  that  the  near-limit-cycle  has  been  achieved.  This  limit- 
cycle  solution  is  an  "ordinary  attractor",  the  attractor  being  the  1-D 
object  to  which  the  phase-space  trajectories  are  attracted  at  all  times. 

The  complete  motion  is  known  once  the  geometry  of  the  attractor  is 
determined.  Hence,  it  is  also  possible  to  compute  the  mean  flow  by 
averaging  the  flow  over  one  complete  cycle;  similarly,  it  is  possible  to 
determine  the  Reynolds  stresses  from  first  principles,  although  these 
computations  have  not  been  performed  in  the  present  study. 

For  the  flow  configuration  with  af-30°,  the  physics  changes 

dramatically,  as  seen  in  Figs.  10a,b,c.  The  curves  for  the  force 
coefficients  show  that  one  limit  cycle  consists  of  two  TE  vortices 


sheddings.  There  are  now  two  shedding  frequencies,  or  modes,  associated 
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with  this  more  complex  attractor.  As  shown  in  Fig.  10a,  the  first  frequency 
is  associated  with  the  shedding  which  takes  place  at  point  1,  whereas  the 
second  frequency  is  associated  with  the  shedding  which  takes  place  at  point 
3.  The  LE  shear  layer  associated  with  the  first  frequency  is  thinner  and 
more  intense,  as  compared  to  that  associated  with  the  second  frequency.  The 
energy  now  oscillates  between  the  two  unstable  modes  through  a  nonlinear 
coupling.  The  appearance  of  subharmonics  signifies  small  modulations  in  the 
shedding  frequency.  This  flow  field,  with  its  two  natural  incommensurate 
frequencies,  is  referred  to  as  a  quasiperiodic  flow,  also  known  as  "Hopf 
bifurcation"  into  an  invariant  torus.  From  Fig.  10a,  it  is  clear  that  the 
initial  state  at  point  10  of  a  new  cycle  is  slightly  different  from  that  at 
point  8.  If  the  phase-space  trajectories  were  drawn,  this  solution  may  very 
well  show  a  tendency  to  fill  a  rather  significant  surface  area  of  the  torus. 
Finally,  it  should  be  noted  that  the  CD  peaks  in  Fig.  10b  correspond  to 

points  2  and  9  in  Fig.  10a. 

For  the  case  with  af-53°t  the  results  obtained  up  to  t-71*  may  be  far 

from  approaching  an  asymptotic  state.  Figure  11a  shows  some  resemblance  of 
quasiperiodic  flow  with  three  incommensurate  frequencies,  this  fact  being 
further  supported  by  the  curves  in  Figs.  11b,  c.  From  their  numerical 
experiments,  Grebogi,  Ott  and  Yorke  [1983]  have  also  shown  the  existence  of 
quasi  period! city  with  three  Incommensurate  frequencies.  The  state  of  the 
system  at  a  given  time  instant  in  one  cycle  is  not  quite  repeated  at  the 
corresponding  time  instant  in  the  subsequent  cycle.  The  phase-space 
portrait,  not  shown  here,  is  very  complex,  where  the  surface  has  folded 
repeatedly  onto  itself,  so  that  it  appears  to  be  a  strange  attractor.  This 
is  an  indication,  although  preliminary,  that  the  flow  may  be  exhibiting  a 
route  to  chaos.  Some  of  the  rigorous  approaches  for  characterizing  a 
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strange  attractor  consist  of  the  determination  of  (i)  the  Lyapunov  exponent; 
(ii)  the  fractal  dimension  of  the  attractor,  which  is  related  to  the  number 
of  degrees  of  freedom;  and  finally,  (iii)  the  Kolmogorov  entropy.  These 
indices  still  need  to  be  studied  thoroughly  in  order  to  rigorously  analyze 
the  route  to  chaos  in  a  meaningful  way. 

The  overall  state  of  the  total  flow  system  was  examined  in  terms  of  the 
frequencies  associated  with  the  various  sheddings  and  an  invited  paper  based 
on  these  results  was  presented  by  K.  Ghia  et  al .  [1985b].  A  written  version 
of  the  paper  was  also  prepared  by  K.  Ghia  et  al.  [1986]. 


Cambered  Joukowski  Airfoil 

The  2-D  unsteady  flow  analysis  was  continued  for  a  cambered  Joukowski 
airfoil.  A  typical  clustered  conformal  C-grid,  with  (230,  46)  points,  is 
shown  in  Fig.  12  for  the  Gottingen  580  airfoil  (i.e.,  an  11.8  percent  thick 
cambered  Joukowski  airfoil  with  a  zero-lift  angle  of  attack  aQ  -  -5.711®), 

at  effective  flow  angle  of  attack  a  -  30®. 

© 

Results  are  obtained  for  three  flow  configurations  with  Re  -  1000  and 
a@  •  5.711®,  15®  and  30®,  the  last  two  being  in  the  stall  and  post-stall 

flow  regimes,  respectively.  Here,  the  effective  angle  of  attack  is 

ae  m  af  ~  ao*  with  being  the  geometric  angle  of  attack  between  the  chord 

and  the  free-stream  direction. 

Figure  13  shows  the  inviscid  starting  solution,  the  grid  distribution 

in  the  near  field  and  the  steady-state  stream-function  and  vorticity 

contours  for  the  case  with  Re  -  1000  and  a  -  5.711®.  The  stream-function 

e 

contours  show  a  mildly  separated  region  in  the  vicinity  of  the  trailing  edge 
(TE).  Laminar  boundary  layers  prevail  on  both  the  suction  and  pressure 
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surfaces,  as  observed  in  the  vorticity  contours,  which  also  show  a  tongue¬ 
like  behavior.  The  streamwise  extent  of  the  separated  flow  near  the  TE  is 
0.15c,  as  compared  to  approximately  0.4c  for  the  symmetric  Joukowski  airfoil 
with  a  -  5°,  as  given  by  K.  Ghia,  Osswald  and  U.  Ghia  [1985b]. 

The  time  history  of  the  lift  and  drag  coefficients  CL  and  CD, 

respectively,  is  depicted  in  Fig.  14  for  the  case  with  Re  -  1000,  o  -  15°. 

© 

A  time  asymptotic  limit-cycle  solution  evolves  by  approximately  t  -  12,  as 
compared  to  t  »  31  for  the  case  of  the  corresponding  symmetric  airfoil. 

K.  Ghia  et  al.  [1985b]  had  defined  the  period  as  the  time  Interval  required 
for  the  Lg-norm  of  the  deviation  of  the  instantaneous  state  of  the  flow  from 

a  reference  initial  state  to  become  smaller  than  a  specified  tolerance.  For 
the  present  configuration,  the  period  is  established  by  examining  the 
successive  maxima  of  CL  in  Fig.  14.  The  period  associated  with  this  limit- 

cycle  (ordinary  attractor)  is  1.046  characteristic  time  units  and  the 

corresponding  Strouhal  number  S  -  fc  (sin  a-)/U  -  0.154;  based  on  a  ,  its 

r  •  e 

value  is  S  -  0.247.  The  corresponding  period  for  the  symmetric  airfoil  is 
e 

1.416  and  S  -  0.18.  Hence,  the  effect  of  camber  is  to  increase  the 
ae 

frequency  of  the  periodic  shedding  of  large-scale  vortices  from  the  TE 
region  and  from  the  separation  zone  on  the  suction  surface. 

The  coefficients  of  lift  and  drag  Cp  for  the  case  with  Re  -  1000, 

a  «  30°  are  shown  in  Fig.  15.  These  curves  show  that,  even  at  t  -  52,  the 
© 

flow  has  not  yet  reached  an  asymptotic  state.  Also,  these  curves  show  that 
one  cycle  consists  of  two  TE  sheddings.  Thus,  there  are  two  shedding 
frequencies  associated  with  this  attactor;  these  correspond  to  the  sheddings 


associated  with  points  1  and  2  in  Fig.  15.  The  LE  shear  layer  associated 
with  the  first  frequency  is  thinner  and  more  intense  as  compared  to  that 
associated  with  the  second  frequency.  This  flow  field,  with  its  two  natural 
incommensurate  frequencies,  is  again  referred  to  as  a  quasiperiodic  flow. 

The  phase-space  portrait  is  complex  and  is  tending  towards  being  a  "strange 
attractor".  The  corresponding  Poincare  sect-^n  indicates  that  the  attractor 
may  be  a  thin  torus.  These  results  are  similar  to  those  of  the  symmetric 
airfoil  (Fig.  11),  except  that  the  history  appears  more  chaotic  and  does 

not  seem  to  be  tending  towards  a  limit  cycle.  Also,  there  is  a  peak  due  to 
the  shedding  from  the  LE  at  point  3;  this  was  not  present  in  the  results  for 
the  symmetric  airfoil.  The  time  instants  at  which  the  detailed  flow  results 
are  presented  correspond  to  the  TE-LE-TE  sheddings. 

The  instantaneous  stream-function  contours,  presented  in  Figs. 
I6a,c,e,g,i,  show  massively  separated  flow,  with  large  eddies  present  over 
the  suction  surface  as  well  as  in  the  wake.  Figures  I6a,g,i  show  the 
presence  of  multiple  separations,  whereas  Fig.  16c  shows  the  presence  of  two 
counterclockwise  co-rotating  bubbles,  aft  of  the  shoulder,  towards  the  TE 
and  these  bubbles  are  in  the  process  of  coalescing.  The  corresponding 
vorticity  contours  are  shown  in  Fig.  I6b,d,f,h,j  and  various  vortex 
interactions  can  be  observed  from  this  figure.  In  Fig.  16b,  the  TE  vortex 
has  Just  been  shed,  a  new  TE  vortex  intensifies  as  shown  in  Fig.  I6d  and  is 
being  Just  separated  from  the  TE  by  the  growing  LE  vortex.  Figure  l6f 
corresponds  to  shedding  of  the  LE  vortex.  The  state  shown  in  Fig.  16J 
corresponds  to  TE  shedding  and  has  features  similar  to  those  in  Fig.  16b. 

The  results  of  this  effort  may  be  summarized  as  follows.  Up  to  the 
stall  regime,  camber  causes  the  flow  fields  to  be  dominated  by  the  TE 
geometry.  Also,  the  extent  of  the  separated  flow  is  diminished  both  in  the 


streamwi.se  and  the  lateral  dimensions  and  the  shedding  frequency  is 


increased.  For  flow  fields  in  the  post-stall  regime,  the  C.  time  history 

L 

becomes  more  chaotic  and  shows  a  peak  associated  with  LE  shedding.  The 
computation  of  the  flow  field  needs  to  be  continued  to  larger  t,  to  predict 
its  further  behaviour  with  definite  assurance.  A  paper  based  on  these 
results  was  presented  by  Osswald  et  al.  [1986]  and  appears  in  the  Conference 
Proceedings. 

Circular  Cylinder 

For  flow  past  airfoils  at  high  angle  of  attack,  the  airfoil  behaves 
like  a  bluff  body.  The  large-scale  structure  of  the  separated  flow  past  a 
bluff  body,  e.g.,  a  circular  cylinder,  is  not  understood  well,  except  at 
extremely  low  values  of  the  Reynolds  nunber  Re  around  120  where  Re  is  based 
on  the  cylinder  diameter.  The  understanding  of  separated  flow  past  the 
model  problem  of  a  circular  cylinder  is  crucial,  as  it  may  lead  to  a  better 
understanding  of  more  complex  flows  such  as  separated  flow  past  airfoils 
with  or  without  solution  bifurcations  associated  with  lift  hysteresis, 
symmetry  breaking,  etc.  Also,  a  thorough  knowledge  of  steady  separated  flow 
could  aid  in  analyzing  unsteady  separated  flow  leading  to  incipient 
transition  and,  eventually,  to  chaos. 

Hence,  the  flow  past  a  circular  cylinder  was  computed,  first  without 
assuming  symmetry.  The  instantaneous  streamlines  and  vorticity  contours  are 
presented  in  Fig.  17  for  Re  -  200.  The  occurrence  of  the  bulge  phenomenon 
for  Re  •  200  and  t  i  1.0,  first  reported  by  Bouard  and  Coutanceau  [1980],  is 
also  computed  satisfactorily.  It  leads  to  an  alteration  in  the  vorticity 
contours  close  to  the  surface  of  the  cylinder  in  the  separated  zone.  For 
the  corresponding  symmetric  configuration  with  Re  -  500,  the  instantaneous 
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stream-function  and  vorticlty  contours  are  delineated  in  Fig.  18  at  t  •  100, 
225,  800  and  1600.  The  eddy  length  L,  as  seen  from  the  stream- function 
contours,  increases  from  approximately  6  diameters  at  t  -  100  to  27 
diameters  at  t-1600;  the  asymptotic  results  of  Fornberg  [19853  show  this 
length  to  be  36  diameters.  Further,  the  region  between  the  cylinder  and  the 
main  eddy  grows  from  about  1  diameter  initially  to  5  1/2  diameters  by  t  • 
1600;  the  corresponding  asymptotic  results  of  Fornberg  [1985]  show  this 
distance  to  be  approximately  8  diameters.  This  region  can  be  considered  as 
consisting  of  two  parts:  the  front  part  where  the  separating  streamline 
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grows  in  width  as  0(Re  ),  followed  by  a  ’transition’  region  in  which  the 
growth  rate  adjusts  such  that  the  width  can  grow  as  0(Re)  in  the  main  eddy. 
Also,  in  this  transition  region,  the  vorticity  increases  from  zero  to  the 
level  of  nearly  uniform  vorticity  present  in  the  main  eddy.  The  kink  in  the 
contours  of  the  vorticity  is  conjectured  to  be  related  to  this  change  in  the 
growth  rate  of  the  width  of  the  eddy.  The  results  for  the  velocity  do  not 
show  0(1)  thickness  for  both  the  shear  layer  past  the  center  of  the  main 
eddy  and  the  return  jet,  as  proposed  by  Smith  [1985];  the  present  results 
need  to  be  scrutinized  further  to  provide  the  accurate  flow  structure.  The 
Navier-Stokes  calculations  show  that,  near  Re  -  500,  the  translating 
cylinder  drags  with  it  a  massive  eddy. 

These  results  far  flow  past  cylinder  were  presented  by  K.  Ghia  et  al. 
[1986]  at  the  I.U.T.A.M.  Confer ence;  a  written  version  of  the  paper  is  to 
appear  in  the  Conference  volume. 

2.4  Unsteady  Three-Dimensional  Navier-Stokes  Analysis 

The  effort  directed  to  this  area  of  study  was  focused  primarily  on  the 
careful  selection  of  the  formulation  of  the  governing  partial  differential 


equations  and  the  proper  choice  of  discretization  to  produce  sparse  matrices 
with  repetitive  block  structure.  This  structure  is  essential  for  efficiency 
of  the  inversion  technique  with  a  high  degree  of  vectorizability  in  the 
solution  phase.  The  goal  was  to  develop  a  direct-solution  technique  for  the 
3-D  unsteady  Navier-Stokes  equations  for  general  geometries.  This  was  based 
on  the  fact  that  direct- solution  methodology  developed  by  the  principal 
investigators  for  2-D  unsteady  Navier-Stokes  equations  has  proved  very 
robust  and  efficient. 

After  much  careful  study,  the  velocity-vorticity  (V,o>)  formulation  has 

been  selected,  as  opposed  to  the  primitive-variable  (V,p)  and  the  several 

vector-potential-vorticity  (A,w)  formulations.  All  vector- potential 
formulations,  which  include  stream-like  functions,  toroidal  and  poloidal 
potentials,  etc.,  suffer  from  the  necessity  for  more  boundary  conditions  on 

the  vector  potential  A  than  are  available  naturally  from  the  physics  of  the 
flow  problem.  This  difficulty  can  be  traced  to  the  fact  that  the  vector 
potential  is  gauge  invariant  and  'extra*  boundary  conditions  are  required  to 
uniquely  select  a  specific  gauge.  Conversely,  the  non-physical  boundary 
conditions  selected  must  be  consistent  with  a  specific  gauge  function; 
otherwise,  numerical  convergence  difficulties  can  be  expected. 

At  first  glance,  the  primitive- variable  (V,p)  formulation  appears  to  be 
computationally  more  efficient,  since  it  requires  the  solution  of  only  four 
unknowns,  V1 ,  V2,  V^,  and  p,  rather  than  the  six  unknowns,  V1 ,  V2,  V^.  u»1 , 

<i>2,  and  uy  occurring  in  the  velocity-vorticity  (V,u>)  formulation.  Indeed, 

much  effort  has  been  directed  by  the  technical  community  at  primitive- 
variable  techniques  for  compressible  flows.  However,  in  the  present  effort, 


it  is  determined  that  the  velocity-vorticity  (V,w)  problem  can  be 
discretized  in  such  a  manner  as  to  produce  a  sparse  matrix  problem  with 
repetitive  block  structure.  As  a  result,  the  computational  work  associated 
with  the  implicit  determination  of  one  unknown  can  be  effectively  eliminated 
from  the  solution  of  the  velocity  problem,  while  the  computational  effort 
for  yet  another  unknown  can  be  eliminated  from  the  vorticity- transport 
equation.  Consequently,  velocity-vorticity  techniques  deserve  attention. 

A  philosophical  point  to  bear  in  mind  is  that  the  (V,w)  formulation 
leads  to  a  more  natural  decoupling  of  the  governing  equations  than  occurs  in 
the  primitive- variable  technique.  Specifically,  pressure  is  an  essential, 
nonlinear  function  of  velocity,  whereas  vorticity  is  a  linear  function  of 

velocity.  Indeed,  the  (V,w)  formulation  can  physically  separate  the  spin 
dynamics  of  a  fluid  particle  (vorticity- transport  problem)  from  the 
translational  kinematics  of  the  fluid  particle  (the  elliptic  velocity 

problem).  This  natural  decoupling,  which  occurs  in  the  (V.to)  formulation, 
may  be  seen  to  directly  translate  into  simple  and  direct  application  of  the 
necessary  boundary  conditions.  Furthermore,  complete  internal  consistency 

can  be  maintained  with  the  (V,u)  formulation,  i.e.,  all  Integral  vorticity 
constraints,  all  Integral  velocity  constraints,  solenoidal  velocity,  etc., 
can  be  algebraically  guaranteed,  irrespective  of  grid  size  and  time-step 
discretization,  throughout  the  entire  flow  solution. 

Vfu  [1985]  and  his  associates,  Fasel  [1976],  and  Gatsky,  Grosch  and 

Rose  [1982]  have  used  the  corresponding  2-D  (V,u>)  formulation  and  are 

presently  developing  3~D  algorithms  using  the  (V,w)  formulation.  The 
present  analysis,  however,  aims  to  employ  efficient  direct  inversion  for  the 


elliptic  velocity  problem  which  can  be  formulated  so  as  to  produoe  a 


uniquely  determined  nonsingular  vector-matrix  problem  with  repetitive  sparse 
block  structire. 

2.5  Pressure-Field  Evaluation  for  Unsteady  Flows 

For  unsteady  flows  analyzed  using  the  («,  formulation  of  the  Navi er¬ 
st  okea  equations,  the  surface  distribution  of  the  pressure  may  be  obtained 
by  simply  integrating  the  tangential  component  of  the  momentum  equation 
along  the  surface.  For  determination  of  pressire  in  the  total  flow  field, 

U.  Ghia  et  al .  [1976]  have  shown  that  a  path- independent  pressure  field 
results  only  from  the  solution  of  the  geumann-Poiaaon  pressure  problem 
formed  by  the  divergence  of  the  aomenttai  equations.  With  Xeiaiann-type 
boundary  conditions,  this  problem  admits  a  solution,  which  is  unique  upto  an 
arbitrary  additive  constant,  only  when  the  boundary  values  and  the  soiree  in 
the  differential  equation  satisfy  Green's  lntep*al  theorem.  This  places 
certain  very  specific  requirements  on  the  discretization  of  the  Neuaann- 
Poisson  problem.  This  work  is  being  pursued  with  the  assistance  of  Collopy 
[1986]  and  comprises  his  Master' s  degree  thesis,  expected  to  be  completed 
shortly.  The  application  considered  was  the  flow  through  an  orifice  in  a 
doubly  infinite  circular  pipe.  The  infinite  extent  of  the  flow  domain  in 
the  streamwlse  direction  necessitated  additional  considerations  in  order  to 
maintain  boundedness  for  the  working  variable  for  the  pressure. 

2.6  Supercomputer  Usage  and  Graphical  Post-Processing  of  Data 

All  of  the  research  pursued  towards  this  y  ant  makes  extremely 
extensive  use  of  computer  resources.  For  acceptable  productivity,  all 
available  computer  facilities  are  sade  use  of,  Judiciously.  During  the 


initial  development  stages,  the  programs  are  tested  and  run  using  the 
facilities  of  the  University  of  Cincinnati  (UC);  these  include  the  Aerospace 
Engineering  Department's  Perkin  Elmer  3250  supermini  computer  system  with 
dual  processors  for  parallel  computation  possibility  and  the  UC  Computer 
Center's  AMDAHL  470  V/7A  mainframe  facility.  Most  of  the  actual  results  are 
eventually  generated  using  the  NASA-Lewis  CRAY  XMP  supercomputer  system  for 
the  IPNS  computations  and  the  NASA-Langley  CYBER  205  supercomputer  system 
for  the  unsteady  Navler-Stokes  computations  (2-D  as  well  as  3-D).  All  four 
of  these  systems  are  equipped  with  their  respective  peripheral  devices  which 
are  employed  for  post-processing  the  results  of  the  computations  and  for 
presenting  these  results  in  suitable  graphical  form.  All  of  the  graphics 
packages  employed  are  written  by  the  present  research  team  and  are  very 
efficient.  However,  the  retrieval  of  the  large  data  bases  from  the  remotely 
located  supercomputers  is  extremely  inefficient  via  the  present  long¬ 
distance  communication  network.  Therefore,  plans  are  presently  underway  to 
purchase,  through  AFOSR  funds,  a  superworkstation  and,  through  the 
University  of  Cincinnati,  have  a  leased  direct-communication  line  to  the 
NASA-Langley  supercomputer  site;  this  is  most  essential  for  the  3-D  unsteady 
flow  simulations,  which  are  presently  being  pursued.  The  limited 
availability  of  supercomputer  time  continues  to  create  a  great  deal  of 
difficulties.  Much  valuable  time  is  often  expended  in  transferring  programs 
from  one  computing  facility  to  another,  in  the  hope  of  expediency,  but  this 
necessitates  continual  changes  in  the  programs  and  ends  up  consuming 
additional  personnel  time.  Easier  availability  of  larger  amounts  of 
supercomputer  time  would  be  most  conducive  for  rapid  progress  in  this 
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FIG.  5.  2-D  COMPUTATIONAL  REGIONS  FOR  APPLICATION  OF  ADI  SCHEME  TO 

HYBRID  GRID  FOR  COMPOSITE  SOLUTION. 
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ENLARGED  HORIZONTAL  SCALE 
■SURFACE  PRESSURE  DISTRIBUTION 
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GEOMETRY  AND 

INFLOW  AND  OUTFLOW  BOUNDARY  CONDITIONS 
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SURFACE  SHEAR  PARAMETER 


FIG.  6  .  RESULTS  FOR  SYMMETRIC  CHANNEL  WITH  EXPONENTIAL  CONSTRICTION; 
R«  -  1500,  t  „  /h  -  0.2,  (201  x  61 )  GRID. 
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Flf.  10.  Inatamanaoua  Aarodynaalc  Coafficianta  at  R«  ■  1.000.  af  •  S0°. 
(a)  Lift  C^,  (b)  Drag  CD,  (c)  Moawnt  C^.  -  Syamatric 

Jotikowald  Airfoil. 
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Fig.  13  GSttingen  580  at  Re  *  1000,  a  »  5.711°;  at  t  »  0:  Inviscid  Stream 
Function  and  Grid  Distribution; at  t  *  18.0  Steady-State  Solution. 
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Fig.  15  Lift  and  Drag  Histories 
G&ttingen  580. 
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INSTANTANEOUS  STREAM- FUNCTION  CONTOURS  VORTICITY  CONTOURS 

nG.  17  DEVELOPMENT  OF  FLOW  PAST  CYLINDER  AT  R«  -  200;  SYMMETRY  NOT  ASSUMED. 
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